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(N '. ABSTRACT 
t-H ■ 

Context. The link between the shaping of bipolar planetary nebulae and the mass ejection activity of their central stars is still poorly 
understood. Appropriately characterizing the evolution of the shells ejected during the late stages of stellar evolution and the interac- 
tion between these shells is fundamental to gain insight into the mechanism of nebular shaping. This must include the study of the 
molecular emission, which tracks the mass-loss history during the late AGB and post-AGB stages, when the nebula is being actively 
shaped. 

Aims. Herschel/HIFI provides an invaluable tool by opening a new window (most of the sub-mm and far infrared range is only acces- 
sible from space) from which to probe warm molecular gas (~50-1000 K). This paper presents a radiative-transfer, spatio-kinematic 
modeling of the molecular envelope of the young planetary nebula NGC 7027 in several high- and low-./ 12 CO and l3 CO transitions 
observed by Herschel/HIFI and the IRAM 30-m radio-telescope, and discusses the structure and dynamics of the molecular envelope. 
Methods. We have developed a code which, used along with the existing SHAPE software (Steff en et al .11201 ih . implements spatio- 
kinematic modeling with accurate non-LTE calculations of line excitation and radiative transfer in molecular species. We have used 
i-C 1 this code to build a relatively simple "russian doll" model to account for the physical and excitation conditions of the molecular 

CU envelope of NGC 7027. 

Results. The model nebula consists of four nested, mildly bipolar shells plus a pair of high- velocity blobs. The innermost shell is 
the thinnest and shows a significant jump in physical conditions (temperature, density, abundance and velocity) with respect to the 
* ■ adjacent shell. This is a clear indication of a shock front in the system, which may have played a role in the shaping of the nebula. 

Each of the high- velocity blobs is divided into two sections with considerably different physical conditions. The striking presence of 
H2O in NGC 7027, a C-rich nebula, is likely due to photo-induced chemistry from the hot central star, although formation of water 
by shocks cannot be ruled out. The computed molecular mass of the nebula is 1.3 M G , compatible with that derived from previous 
works. 
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. 1. Introduction less to accurately probe molecular gas with temperatures over 
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100 K. 



^r^ The physical mechanisms behind the shaping of proto-planetary 

CN ' nebulae (PPNe) and vounsnlanetarv nebulae (PNe) has been and Recently, as part of the guaranteed-time key program 
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nebulae (PPNe) and young planetary nebulae (PNe) has been an d 

still is a subject of intense debate (see e.g. lBalick & Frankl2002l) . HIFISTARS, a sample of ten PPNe and young_PNe were ob- 

In particular, PPNe and young PNe usually show high-velocity served wlth Herschel/HIFI d Buiarrabal et al. | |2012|). Although 

outflows, usually in the polar directions, along with slower com- HIFISTARS only includes observations towards the center of 

ponents whose velocities are comparable to those of the winds these nebulae, the kinematics and the excitation conditions of 

£3 1 of AGB envelopes. These fast-bipolar outflows are thought to be the warm molecular gas of the target can be studied with un- 

" the result of the propagation of shocks in the polar regions of the precedented detail thanks to the high spectral resolution of the 

dense AGB shell during the evolution towards PPNe. However, profiles. 

many aspects of this process are poorly understood at present. This work focuses on one of the targets observed within 

Tracking the mass ejecta of planetary nebulae (PNe) and the HIFISTARS program, NGC 7027 (PN G084.9-03.4). With 

proto-planetary nebulae (PPNe) by means of molecular emis- a kinematical age of just 600 years (b ased on its radio flu x, see 

sion lines can provide crucial clues for better understand- lMassonlfl989h and located at 1 kpc (IZiilstra et al]l2008l) . this 

ing this topic. Herschel/HIFI provides an invaluable tool, as compact and young PNe is one of the brightest nebulae in the 

it probes warm molecular gas (~50-1000 K) by producing sky and the most extensively studied PNe so far. NGC 7027 is 

high-resolution spectra in high-excitation transitions (e.g. 12 CO a carbon-rich nebula with a very high-excitation spectrum (e.g. 

7=6-5, 10-9 and 16-15), some of which are unobservable from it shows ion lines such as O iv and Mg v). It hosts one of the 

the ground, rendering ground-based telescopes esentially use- hottest central star s known to date, with a T eff ~200,000 K (e.g 

lLatter et all l2000l) . A small, essentially ellipsoidal, e xpanding 

* Based on observations made with the Herschel Space Telescope ionized shell surrounds the central star dMassonl l 1989b. Further 

operated by ESA and the 30-m radiotelescope operated by IRAM. outwards, a thin H2 shell indicates the presence of a photo- 
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dissociation region (PDR) an d shows signs o f recent interac- 
tion with collimated outflows (ICox et al.ll2.002h . The nebula also 
shows molecular emission from CO beyond the PDR, extending 
further out to 15-20 ar csec from the central star. In particular, 
iNakashima et al.l d2010l) carried out interferometric observations 
in the CO 7=2-1 and 1-0 lines and found that the CO envelope 
could be roughly modeled as an ellipsoidal shell, despite some 
deviations caused by a high-velocity component. The kinemat- 
ics of this cold, molecular envelope seems to be closely linked 
with that of the in ner ionized shell and the PDR (ICox et al.lll997l 
lLatter et alfeoOOl) . 

The amount of molecular mass inferred from the CO emis- 
sion is very significant. In fact, despite hosting one of the hottest 
central stars known, the ionized mass only represents 2% of the 
total measured mass of the nebula, 1 .4 M , while the vast major- 
ity of the nebula (85%) consists of molecular gas, the remainin g 
13% being in the form of neutral atomic gas dFong et al1l2001l) . 
It seems clear that any approach towards understanding the shap- 
ing mechanism of this nebula must take into account the emis- 
sion from molecular gas. 

In this work, we present a detailed study of the molecular en- 
velope of NGC 7027, combining observational data from low-/ 
transitions, observed with the IRAM 30-m radio-telescope with 
data from high-/, observed as part of the HIFISTARS program. 
We hav e developed a code , shapemol, which, used along with 
SHAPE dSteffen et al.ll201 lb . implements spatio-kinematic mod- 
eling with accurate calculations of non-LTE line excitations and 
radiative transfer in molecular species. The high quality of the 
HIFI data, together with shapemol, have allowed us to study, 
for the first time in such level of detail, the kinematics and ex- 
citation conditions of the warm, molecular (CO-rich) gas of a 
PN. 

A preliminar ver s ion o f this work was reported in 
ISantander-Garcia et all d201 ll) . 



2. Observations 

The observations used in this paper (see Fig. [TJ for the modeling 
of the mo lecular envelope o f NGC 7027 have been already pre- 
sented bv IBuiarrabal et all d200lh. 12 C O and 13 CO /=l-0 and 
2-1 data, and IBuiarrabal et al.l (1201 2h . rest of the data consist- 
ing of the high-/ transitions. Here we will briefly summarize 
the main observational characteristics of the data (see Table 1), 
as they have been described in full detail in these papers. CO 
/=l-0 and 2-1 spectra were obtained in 1999 and 2000 us- 
ing the IRAM 30-m at Pico de Veleta (Spain). This is a very 
well characterized ground-based instrument, in which the main 
source of uncertainty in the absolute calibration of the data 
comes from the correction for the opacity of the atmosphere. 
Although this effect is measured every 20 to 30 min, in standard 
conditions we cannot guarantee accuracies better than 20-25% 
when comparing data from Pico de Veleta with those from other 
telescopes, as this is the best accuracy our group has found on 
long-term observations of the same standard sources, specially 
for observations dating a decade or more. Spectra at higher fre- 
quencies are taken from the Herschel/HIFI HIFISTARS GTKP 
(P.I. V. Bujarrabal). In this case, the main source or uncertainty in 
the absolute calibration of the data of this space-based telescope 
comes from the fact that the gain ratio between the two sidebands 
of the receivers is not well known. For the HIFISTARS data we 
have adopted uncertainty values from 15% to 30% depending on 
the receiver band, following the description of the HIFI instru- 
ment c haracteristics and performances given bv lRoelfsema et all 
d2012l) . In all cases we have used the / m b intensity scale, that can 



be directly compared with the model results when these take into 
account the shape of the beam of the telescope (see below). 

In all cases, the spectra are taken towards the central star, lo- 
cated at J2000 R.A. 21 h 07 m 01 s .59 Dec.+42°14'10:'2. The neb- 
ula is sometimes larger than the beam of the telescope (12-22 
arcsec for the 30-m, and 12-33 arcsec for HIFI), and there- 
fore by observing just the central position we do not observe 
its full emission. This is particularly true for the higher frequen- 
cies in both telescopes (at 220-230 GHz for IRAM 30-m and 
1800 GHz for Herschel/HIFI), but this effect has been properly 
taken into account in our modeling (see section 4). For this pur- 
pose we have assumed that the telescope beam shape is well 
represented by a symmetrical 2D-Gaussian solely characterized 
by its FWHP, what has been taken to be equal to the measured 
HPBW of the two instruments at the corresponding frequencies 
(see Table 1). 

The spectral resolution of the instruments is also limited. For 
the IRAM 30-m observations, the data were taken from the (now 
decommissioned) 100 kHz filter-bank, after degrading the reso- 
lution to 300 kHz by averaging every three channels, resulting in 
velocity resolutions of 0.8-0.9 km s for the 1-0 lines, and of 
0.4 km s for the 2-1 lines. For the Herschel/HIFI data, we have 
used the Wide Band Spectrometer (WBS) in standard configura- 
tion, providing a spectral resolution of about 1 . 1 MHz, slightly 
varying across the observed band. The data have been oversam- 
pled in a homogeneous 0.5 MHz spectral channels, that we have 
degraded down to resolutions equivalent to 0.4 to 1.6 km s~') 
depending of the strength of the detected line and the required 
S/N ratio. 

3. Modeling 

The main goal of our modeling is to determine the physical con- 
ditions (densities, temperatures and abundances) and the general 
kinematic structure of the nebula. Note that an accurate determi- 
nation of the 3-D structure from our data is prevented by the lack 
of spatial resolution. In fact, the only information about the spa- 
tial structure is concealed in the frequency-dependent beam size, 
which weighs the emission from the central region and that com- 
ing from the outer region in a different fashion for each transi- 
tion, and which is, for the highest-/ transitions, smaller than the 
molecular CO shell at lower-/ described by previous works (see 
beam size column in Table [1]). Hence, the general geometry of 
the nebula has been deduced both from the analysis of previous 
works and a series of reasonable deductions from the Herschel 
and IRAM data-set (see section |4~TT ). In any case, the following 
model should be considered as a first, simplified approach to a 
much more complex reality. 

The modeli ng process is similar to standar d spatio-kinematic 
modeling (e.g. ISantander-Garcfa et al.ll2004b . except that it in- 
corporates full radiative-transfer solving for different transitions, 
creating synthetic profiles to be matched against the observa- 
tions. We have achieved this by using a sp e cial ve rsion (4.5 1/3) 
of the SHAP£] software by ISteffen et all d2011l) , specifically 
customized for usage alongside our own GDL/IDL-based code, 
shapemolQ. 

SHAPE is a software tool for building spatio-kinematic mod- 
els of nebulae. It allows to easily implement a 3-D structure 
and 3-D velocity field to describe the model nebula, and gen- 
erate synthetic images, position-velocity diagrams and/or chan- 
nel maps for direct comparison with observations. Its versatility 



SHAPE is available on http://bufadora.astrosen.unam.mx/shape/ 
If interested in using shapemol, contact the first author. 
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Table 1. Observing log and profile features. 



Species 


Transition 


Frequency 


Resolution 


T mb peak 1 


r.m.s. 


Integ. flux* 


HPBW 






(GHz) 


(kms 1 ) 


(K) 


(K) 


(Kkms 1 ) 


(") 


IRAM 30 


-m observations 11 














l2 CO 


7=1-0 


115.271 


0.8 


11.9 


7.7xl0~ 2 


332 


22 


12 co 


7=2-1 


230.538 


0.4 


30.9 


9.9xl0~ 2 


667 


12 


l3 CO 


7=1-0 


110.201 


0.9 


0.3 


9.4x1 0" 3 


7.9 


22 


l3 CO 


7=2-1 


220.399 


0.4 


1.3 


3.0xl0~ 2 


31 


12 


Herschel/HIFI observations' 3 




12 co 


7=6-5 


691.473 


0.4 


5.1 


1.3xl0- 2 


137 


31 


12 co 


7=10-9 


1151.985 


0.4 


5.5 


5.5xl0~ 2 


146 


20 


l2 CO 


7=16-15 


1841.346 


0.5 


2.5 


9.3xl0- 2 


52 


12 


l3 co 


7=6-5 


661.067 


0.5 


0.3 


8.6xl0~ 3 


8.3 


31 


l3 CO 


7=10-9 


1101.350 


0.5 


0.2 


2.5xl0- 2 


5.7 


20 


c 18 o 


7=6-5 


658.506 


0.9 


0.03 


6.1X10- 3 


0.7 


33 


OH 


2 n, /2 7=3/2--l/2+ 


1834.75 


1.6 


0.16 


4.8xl0~ 2 


6.3 


12 


p-H 2 


7r a ,j: c =li,i - .o 


1113.343 


0.9 


0.1 


1.9xl0" 2 


2.1 


20 


o-H 2 


7^ f( ,A' f - = ll.O - lo.i 


556.936 


0.5 


0.1 


5.5xlO- 3 


2.1 


39 



Notes: 

1 Maximum intensity of the profile in main beam temperature 7" m b(K), defined as the peak intensity minus the r.m.s. in the spectrum 
* Integrated flux of the profile. 
a Observations published in 1 
b Observations published in 1 



Buiarrabal et al 


2001 


Buiarrabal et al 


2012 



has made it the standard tool in the field, where it has proven to 
be valuable for producing accura te spatio-kinematic descriptions 
of many planetary nebulae (e.g. ISteffen et alj|201 ll iJones et ail 
l20ldlVelazquez et al.ll201~uV 

Although SHAPE implements radiative-transfer solving for 
atomic species, in which the absorption and emission coeffi- 
cients are easily predictable over a large range of temperatures, 
densities and abundances, it is unable, on its own, to work with 
molecular species, either in the thermalized or non-thermalized 
cases. The main reason behind this is the strong dependence of 
the absorption and emission coefficients of each grid cell, not 
only on the local abundance, temperature and density, but also 
on the size of the emitting structure and the velocity of the grid 
cell along the line of sight. 

We designed shapemol to fill this gap. In the following sub- 
section, we describe this software and its role in the modeling 
process. 

3.1. shapemol 

shapemol is a GDL/IDL-based code built to be used as a com- 
plement for SHAPE, allowing for radiative transfer solving in 
molecular transitions. In particular, it performs non-LTE calcula- 
tions of line excitations (based on the well-known LVG approx- 
imation, see below) to compute the absorption (k v ) and emission 
(J v ) coefficients of each individual cell in the grid, for the desired 
transition and species. 

In spectral lines, the values of k v and j v are given by the pop- 
ulations of the involved energy levels and the line profile shape. 
In our case (low-frequency transitions requiring relatively low 
excitation), the profile shape is given by the local velocity dis- 
persion due to thermal dispersion or turbulence. The level pop- 
ulations depend on the collisional transition rates and the radia- 
tive excitation and de-excitation rates, which in turn depend on 
the amount of radiation reaching the nebula point we are con- 
sidering at the frequency of the line (averaged over angle and 
frequency within the local line profile). The calculation of such 
averaged radiation intensity requires a previous knowledge of 
the absorption and emission coefficients in the whole cloud, in 
order to solve the radiative transfer equation in all directions and 
frequencies. Of course, the populations of a high number of lev- 
els must be calculated simultaneously, since the population of 



each one depends on those of the others, and in all the points 
of the cloud. So, the solution of the system becomes extremely 
complex in the general case. 

The problem is greatly simplified when there is a Large 
Velocity Gradient (LVG) in the cloud, introducing important 
Doppler shifts between points that are sufficiently far away. 
When this shift is larger than the local velocity dispersions, the 
points cannot radiatively interact at large scales, so the radiative 
transfer is basically a local phenomenon. The excitation equa- 
tions can be then solved, and k v and j v calculated in each point, 
independently of the rest of the cloud. In any case, the level pop- 
ulations depend in a complex way on the (local) physical con- 
ditions, and the solution requires an iterative process. A very 
unders t andab le version of the general formalism can be seen in 
ICastorl (1 1 970h . The LVG approximation includes the main ingre- 
dients of the problem (collisional and radiative rates, trapping 
when opacities are high, population transfer between different 
levels, etc.) and gives fast, sensible solutions. These excitation 
calculations are quite accurate, even when the LVG conditions 
are barely satisfied, at least for the case of molecular lines from 
shells around evolved stars (e.g. Bujarraba l et al. in preparation; 
iTevssier et al.ll2006l iRamstedt et al.ll2008l) . The approximation 
itself is not necessary to derive the resulting line profiles, which 
can be calculated solving the standard transfer equation using 
the level populations derived from the LVG approximation and a 
local velocity dispersion, as indeed we have done using SHAPE. 

According to the LVG approximation, k v and j v depend in a 
heavily non-linear way on the density n and the temperature T, 
and almost linearly on the product £X, where r is the distance 
of a given point (or code cell) to the central star, V its velocity 
and X the abundance of the desired species. Besides these three 
parameters, the LVG results depend on the logarithmic velocity 
gradient, e = dV/drr/V, but only slightly. We will assume e = 1, 
i.e. a linear dependence of V on r. In this case, the escape proba- 
bility of a photon is given by a simple analytical function of the 
opacity and its calculation is considerably simplified. Such ve- 
locity fields are found in many young planetary nebulae, which 
are basically expanding at high velocity following a 'Hubble' 
velocity law. 

The approach of shapemol consists of relying on a set 
of pre-generated tables of k v and j v as functions of n and T, 
each table corresponding to a value of the product. The 
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values of k v and j Y in each table are individually computed 
for each pair of n and T by iteratively solving a set of equa- 
tions via the convergence algorithm typical of LVG codes. In 
this calculations 20-25 rotational levels were taken into ac- 
count, allowing accurate calculations for temperatures up to at 
least 1000 K, and we used collisional rates fro m the Leiden 
Atom ic and Molecular Database (LAMDA; see ISchoier et alJ 
I2005L http://www.strw.leidenuniv.nl/moldata/). For a set of phys- 
ical conditions, shapemol selects the table with the closest value 
of ^X. Once a table has been selected, shapemol computes k v 
and j y by linear interpolation between the values for the two 
adjacent tabulated values of n and T. The steps in n and T are 
small enough to guarantee that a 1 st order interpolation is a good 
approximation between two consecutive values. Finally, given 
the roughly linear dependence of k v and j v on the product ^X, 
the software scales the computed absorption and emission co- 
efficients according to the ratio of the desired value of ^X to 
that of the selected table; to avoid significant errors, calculations 
were performed for a large number of values of this parameter. 
As mentioned, the absorption and emission coefficients calcu- 
lated in this way for each cell of the model nebula are used to 
solve the standard radiative transfer equation and calculate the 
line profile leaving the nebula in a given direction; at this level, 
we must assume values of the local velocity dispersion in each 
point. 

The current version of shapemol is available from the first 
author. At its present state, it is able to compute the 7= 1-0, 
2-1, 6-5, 10-9 and 16-15 transitions of the 12 CO and 13 CO 
species. The values of n and T in the available pre-generated ta- 
bles range from 10 to 1000 K in steps of 5 K for the temperature, 
and from 10 2 cm~ 3 to 10 7 cm 4 in multiplicative factors of vTo 
for the density. A more refined version, which will include more 
species and transitions, is under progress and will be presented 
by Santander-Garcfa et al. (in preparation). 

3.1 .1 . Modeling process using shapemol 

The modeling process using SHAPE and shapemol is as follows: 

1. We model a nebula with several distinct outflows or struc- 
tures using the standard tools in SHAPE. The nebula is then 
diced in a large number of cells (~tens of thousands) using a 
3-D grid. Each cell in this 3-D grid is characterized by a iden- 
tification code (unique per each structure), a position, a ve- 
locity, a density and a temperature. These values are dumped 
into an ascii text file, one line per grid cell. 

2. shapemol reads the text file and let us select the desired 
molecular species and transition, the molecular abundance of 
each different structure and a value for the micro-turbulence, 
<5v- The software then computes the absorption and emission 
coefficients using the LVG approximation and adds them to 
the text file along with the frequency of the selected transi- 
tion and the micro-turbulence value. 

3. SHAPE reads the text file generated by shapemol and per- 
forms full radiative transfer solving for the whole 3-D grid, 
creating a data-cube consisting of a stack of images of the 
nebula, one per frequency resolution element (i.e. spectral 
channel). The radiative transfer solving used in this step is 
exact and implements no approximations. This data-cube 
is then convolved, channel by channel, with the telescope 
beam, which is simulated by a Gaussian with a FWHM equal 
to the telescope beam HPBW. The profiles from the cen- 
tral four pixels are added, and the intensity / in each chan- 
nel is divided by the projected area of the four cells on the 



plane of the sky. A final profile is then shown, in units of 
W s Hz 1 rrT 2 str -1 , which once converted to T m \, allow 
for direct comparison with observations. 

3.1.2. Numerical errors using SHAPE+ shapemol 

We tested SHAPE+shapemol by testing its results against ana- 
lytic theoretical predictions in a spherical nebula under several 
circumstances of abundances, thicknesses, turbulence, velocity 
distributions and beam sizes. In particular, we modeled a sta- 
tionary nebula with an abundance range large enough to make it 
optically thin or thick; a nebula expanding at a constant veloc- 
ity and another one following a Hubble-like expansion pattern, 
for which analytical predictions can be made with the help of 
the LVG approximation. Also, we explored the flux measured 
off-center with a small beam size in all the previous scenarios 
and compared it with the theoretical predictions. All the previ- 
ous tests were performed for different values of the abundance, 
velocities and micro-turbulence values, always fulfilling the re- 
lation V > <5y, as the LVG approximation demands. 

In all cases, the error of the model with respect to the pre- 
dicted intensity was <10%. This error includes numerical errors 
due to the grid size, limited by the RAM memory of the system, 
as well as the rounding errors of computational algorithms. 

A more detailed discussion of the errors and the performed 
tests will be given in Santander-Garcfa et al. (in preparation). 

4. Results 

In this section we discuss the general structure and dynamics we 
assumed as premises to our modeling (section l4~TI ). and describe 
in detail the morphology, kinematics and excitation conditions of 
the resulting best-fit model for the emission in 12 CO and 13 CO 
(section l4~2l . as well as its uncertainties (section l4~3l l. Finally, we 
provide a qualitative analysis of the emission in other molecules 
(section l4.4t . 

4.1. Adopted general structure and dynamics 

The general structure and dynamics of the m odel nebula have 
been deduced from previous imaging (e.g. iNakashima et al.l 
2010, Cox et al. 2002) and several general properties of our data. 
INakashima et al.l (120101) find that the molecular gas is located in 
a hollow axisymmetric shell in expansion, reaching out to 8-10 
arcsec from the central star in the equatorial directi on. This is 
slightl y more extended than the H2 shell imaged by ICox et all 
d200l . which reaches out to ~5 arcsec in the same direction. 
The molecular shell is also known to sho w a low-brightness rel- 
atively extended halo (iFong et al.ll2006l) . reaching out to 20-23 
arcsec from the central star in the equatorial direction. 

Fig. Q] shows the observational profiles of the 12 CO 7=1-0, 
2-1, 6-5, 10-9, 16-15, and 13 CO 7=1-0, 2-1, 6-5 and 10-9 
emission (black histogram), along with the synthetic profiles of 
the best-fit model (red line, see section 14/21 . There are a num- 
ber of features we can deduce from these observational profiles 
to set the basis for our model. In particular, the blueward in- 
tensity fall of the 12 CO 7=6-5 transition is indicative of strong 
self-absorption caused by high-opacity and stratification in tem- 
perature, which would decrease outwards; also, the peaks of the 
emission at low and high 7 occur at different velocities, with 
absolute radial velocity (with respect to the systemic velocity) 
increasing with 7, which in turn indicates simultaneous stratifi- 
cation in velocity and temperature. 
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We have associated the highest-velocity, low-intensity fea- 
tures visible in our data on both sides of the center of each pro- 
file with a pair of unresolved features visible at those velocities 
in the 12 CO 7=2-1 maps by iNakashima et alJ d201dh . They are 
confined in two small regions along the projection of the nebular 
main axis on the plane of the sky, at offsets of ~4 arcsec. They 
appear to lie inside the main envelope of the nebula, closer to the 
central star. We have interpreted these as a pair of high-velocity 
blobs ejected by the central star in the polar direction. 

Some features and parameters of the model can be deduced 
exclusively from the conclusions of previous works. Given that 
the profiles from Herschel/HIFI and the IRAM 30-m radio- 
telescope contain little information on the morphology of the 
molecular envelope, we have assumed a ba sic geometry and size 
simila r to those found bv lCox et al.l d2002l) and lNakashima et al] 
(l2010h for H2 and CO respectively, consisting of a bipolar 8- 
shaped shell with a wide waist. We have assumed the inclination 
to the line of sight and position angle of the main axis of ou r 
model to be the same as those found by Na kashima et all (l2010h . 
i.e. /=60° and P.A.=155 °. Summarizing , we have extended the 



structure d escribed bvlCox et al.l (120021) taking into account the 
findings of INakashima et alf ( 2010lK 

The co mparison betwee n the position-velocity plots in H2 
(Fig. 4 in ICox et all 120021) and in 12 CO 7=2-1 (Fig. 5 in 
INakashima et al J2010h . confirms the stratification of the temper- 
ature and velocity. The velocity width of the IT envelope, which 
lies closer to the star than the CO envelope, is ~15 km s _1 larger 
than that of 12 CO 7=2-1, making it reasonable to assume that 
the innermost shell of the CO envelope is also the hottest, and 
that the shells get cooler and more diffuse with increasing dis- 
tance to the star. 

In our modeling, we have adopted a distance of 1 kpc follow- 
ing the results obtained by Zijlstra (2008) from the comparison 
of the velocity field and apparent expansion at radio wavelengths 
during a 25-year monitoring period (980+100 pc). 

We consider the aforementioned general considerations and 
parameters as fixed in our modeling. Additionally, for the den- 
sit y, tempera t ure, a nd abundance, we used the values found 
by ICox et al ] (I2002h in their Fig. 13 as starting points for the 
model iteration process. Note, however, that our model focus 
on studying the general physical properties and will necessar- 
ily be simple compared with the complex level of details (e.g. 
holes, clumping, inhom oge neities, etc .) show n by the data in 
INakashima et all d2010t) and lCox et all d2002l) . The model will 
therefore focus on reproducing the main features of the observa- 
tional data, but will not provide an accurate fit of the inhomo- 
geneities and other local details. 

4.2. Best-fit model 

The model was fit to the data following an iterative trial-and- 
error procedure. Starting from a set of initial assumptions based 
on previous works (see previous subsection), we modified each 
parameter over a large range and solved the radiative transfer 
equation to obtain a synthetic profile for each of the nine 12 CO 
and 13 CO observational transitions. Each of these synthetic pro- 
files were overlaid on the corresponding observational profiles, 
and the global quality of the fit assessed by eye as a compro- 
mise between the shape of the profile and the intensity at each 
frequency. When a fair fit was found for all the transitions, we 
varied the parameters of the model over a smaller range which 
still provided a fair fit, in order to find the uncertainty of each 
parameter. This process was repeated several times starting from 



a different set of parameters. In all cases, the model converged 
to similar values of every parameter. 

A schematic representation of the best-fit model can be seen 
in Fig.|2l while the associated parameters are shown in tables [2] 
and [3] Note that due to the vast difference in transition frequen- 
cies and origin of the data (ground and space-based), we cannot 
rule out the presence of relative calibration errors. Therefore, we 
have included a intensity factor per transition, i.e. the factor to 
be applied to the intensity of the modeled profile. This factor 
also accounts for computational errors, and is to be considered 
as the conservative error of the intensity of our model. In all 
cases, the intensity of the model profiles were within a 25% of 
the observed profiles (see the applied intensity factors in Fig.[TJ. 
It is noteworthy to point out that the predicted peak intensity of 
the model at the 13 CO 7=16-15 transition, 30 mK, is compatible 
with the non-detection of such line, which peaks at 50±30 mK, 
at a resolution in velocity of 4 km s . 

In the following we give a detailed description of the differ- 
ent features of the model, i.e the main body of the nebula and the 
high-velocity polar blobs. These results are essentially compat- 
ib le with, but more detailed tha n the preliminary ones provided 
in lSantander-Garcia et alj|201 ll whose analysis did not include 
13 CO transitions. 

In all cases, we found the best-fitting values of the 
12 CO/ 13 CO abundance ratio and micro-turbulence 6y to be 50, 
and 2 km s -1 , respectively. The best-fit for the systemic velocity 
was Vlsr=26.3 km s . 

4.2.1 . Main body of the nebula 

The main body of the nebula consists of four nested, wide waist, 
8-shaped shells, arranged in a fashion similar to a Matryoshka 
or Russian Doll, with each shell in contact with the adjacent 
ones. Their s hape is roughly s imilar but more extended than 
the model by ICox et ail d2002l) for the H2 emission. Each of 
the shells is characterized by a relatively simple set of param- 
eters: equatorial and polar maximum radii, r e and r p , a thick- 
ness Ar, single values of the abundance X(CO), temperature 7, 
and density n for the whole shell; and two characteristic ve- 
locities, V r and V z . The total expansion velocity is the combi- 
nation of these velocities, where V r is constant, outwards from 
the central star in the radial direction, and V z is constant, di- 
rected along the polar axis, also outwards from the equator. 
In all cases, V z is only triggered for distances to the equa- 
tor larger than 1.5xl0 16 cm. A similar behavior of the veloc- 
ity is observed in othe r young PNe or PPNe such as M 1-92 
dBuiarrabal et aljl998l). M 2-56 dCastro-Carrizo et al.l2002l) and 



CRL 618 dSanchez Contreras et al .1 12004 



Each shell is tagged by an index according to their location: 
I, Ml, M2, O, after inner, middle 1, middle 2 and outer, respec- 
tively. The inner shell I is located ju st outside the PDR , surround- 
ing the H2 emission shell model bv lCox et alJ d2002l) . and is the 
main contributor to the features seen in the highest-7 spectrum, 
while the outermost shell O is the coldest, extending to a max- 
imum distance of 2.25xl0 17 cm, which projects to 15 arcsec 
on the plane of the sky. This is sligh tly larger than the struc - 
ture visible in 7=2-1 in the maps by INakashima et aTl d2010l) . 
bu t similar in size to the envelope revealed by the 7=1-0 maps 
bv lFong et al.l d2006h . This discrepancy can be explained by the 
large extent and smoothness of the envelope, which would make 
this outermost shell prone to flux-loss (as it would be partially 
resolved out) in interferometric observa tions, specially at higher 
frequencies, as INakashima et al.l d2010t) expects for the 7=2- 1 



emission. 
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Fig. 3. The different parameters of the model of the molecular 
envelope of NGC 7027 against the distance to the central star, for 
the equatorial (thick green line) and polar (dotted-dashed blue 
line) directions, and the polar blobs (dotted red line). The error 
bars have been drawn at the center of each shell for the blobs and 
polar direction only (the errors in the equatorial direction are the 
same as in the polar direction). Note: The velocity (and its error 
bars) of the polar blobs has been divided by 5 for displaying 
purposes. 



Fig. [3] shows the behavior of the different parameters along 
the equatorial and polar directions, while the right side of Fig. 
|2] displays a vector plot of the velocity field. Aside from the 
highly direction-dependent velocities, and given that the shells 
are characterized by single values of the parameters, the behav- 
ior of X, T, and n are similar in every direction, as they experi- 
ment the same jumps, only at different distances from the central 
star. The different parameters are essentially compatible with the 



smooth gradients which would be expected in post-AGB e jecta 
(see e.g. lJusttanont et al.|[l994llGoldreich & Scovillell 19761) . ex- 
cept for the steep jump between the thin, I shell and the thicker 
Ml shell. The abundance increases by a factor 1.25 across the 
I to Ml shell discontinuity, while the density and temperature 
drop by factors 3.75 and 4 respectively. The velocities V r and V z 
both drop by 25%, a change much more abrupt than the slight, 
steady increase in velocity from Ml to the outer layers M2 and 
O. 

We take these jumps as clear indications of a shock front 
traveling outwards through the nebula, in a similar way as the 
low-velocity shock fro nt found for instance in CRL 618 (e.g. 
iBujarrabal et al.l 120101) . in which a double shell is expanding 
against the nebular halo. The post-shock gas would then natu- 
rally be denser, hotter, faster, and less rich in molecules, as a 
small fraction of those would have been dissociated by the phys- 
ical conditions at the outer boundary of the PDR and the pass- 
ing shock front. The presence of this shock implies an increased 
chance of molecular reactions, which could partly explain the 
presence of H2O (see section 4.4). Also, the sh ock would have 
strong implications in the shaping of the nebula dBalick & Frankl 
120021 

While this relatively simple model does not in- 
clude a detailed description of the geometry or inh omo- 
genei ties/peculiarities of the shells (Nakash ima et aljfeOlOl) . it is 
useful to estimate the amount of matter that emits under certain 
physical conditions. Taking into account the volume and density 
of each shell, we can derive the molecular mass distribution. 
Adding up every shell, the total molecular mass we estimate for 
the main molecular envelope of NGC 7027 is 1.15 M n , which 
agree s well with previous estimates of 1.2 M Q by iFong et al.l 
|200TI 

The kinematic age of the inner, presumably shocked shell is 
~ 1700-2000 years (computed at the equator and poles respec- 
tively). The rest of the envelope ranges from ~2800-3000 years 
for the Ml shell to -4300-4700 years for the O one. This time- 
difference is not model-dependent, as it arises directly from the 
obse rvations, and is also visible in the data fr om previous works 
(e.g. lNakashima et al.l2010tlFong et al.l2006l) . This suggests that 
the expansion pattern of the molecular envelope is not compat- 
ible with a single-event, pure ballistic expansion or Hubble-like 
flow, and was instead ejected during a time-interval larger than 
1000 years. It is interesting to note that, in this case, the value 
of n x r 2 (density corrected for the dilution into the medium) de- 
creases with the kinetic age of the outflow, while the expansion 
velocity increases. This is expected if along the AGB evolution 
the stars approaches increases in its mass loss. Another possi- 
bility, though, would be that the nebula has been affected by a 
series of front shocks, resulting in different apparent kinematic 
ages. 

4.2.2. High-velocity polar blobs 

Unfortunately, data from previous works lack enough spatial res- 
olution to resolve the geometry of the polar blobs. Therefore, 
we chose to model them in a simple way as two small cylin- 
ders being ejected from the central star along the nebular main 
axis. The reason for choosing cylinders and not other geometri- 
cal shapes (e.g. spheres) is for mere convenience: each cylinder 
is easily split in two sections with different physical conditions. 
Each section is characterized by a height, a radius, a distance 
from the equatorial plane to its center, and single values of T, X 
and «. The velocity of each section is constant and its direction 
is axial (i.e. parallel to the main axis of the nebula). The best- 
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Table 2. Best-fit model parameters for the main molecular envelope of NGC 7027. The sizes correspond to a distance to the nebula 
of 1 kpc. The uncertainty range is included for every parameter except the sizes of the shells, which depend on the corresponding 
uncertainties on the shell thicknesses. 



Shell 


polar 

(10 16 cm) 


r 1 

equatorial 

(10 16 cm) 


Ar* 
(10 16 cm) 


V? 
(km s" 1 ) 


v; 

(km s 1 ) 


X( n CO) 


n 

(cm" 3 ) 


T 

(K) 


I: Inner shell 
Ml: Middle 1 shell 
M2: Middle 2 shell 
O: Outer shell 


13.5 
15.75 
18.0 
22.5 


7.5 
9.75 
12.0 
16.5 


Q 75+ 0.15 
o. /j_ 015 

9 9 c +0.25 

-0.35 

9 9^+0.35 

^•^ -0.35 

A 5+0-6 


13 +0 - 3 

"■"-0.3 
Q 7^+0.35 
y - -0.25 

io.o:°j 

i i Q+0.4 


7 c+0.4 
'• J -0.5 

5 6 +03 

6 +0 - 4 

7 +0 - 5 

' -"-0.4 


s :::xio 5 
l+^xlO- 4 
41^xl0- 4 
6.7^xl0- 4 


1.5![Jfxl0 5 
4.0^|xl0 4 
2.5+°gxl0 4 
5.0!^xl0 3 


400+ 40 , 
100^° 
601 

25 + J 2 



Notes: The X( u CO) abundances are 50 times lower than their X{ l2 CO) counterparts. 
^ Distance from the star which traces the shell's outer edge, both at the pole and equator. 

* Thickness of the shell. 

* Expansion velocity radially outwards from the central star. 

* Additional velocity component along main axis. Activated for distances to the equator greater than 1.5X10 16 cm. 



fit parameters of the modeling are shown in Tableland Fig. [3] 
(dotted red line). 

The aft section, located closer to the central star, is thin- 
ner, hotter and denser than the fore one, which center is 
7.65xl0 16 cm offset from the star, so that its projection unto the 
plane of the sky coincides with the small, high-veloc i ty reg ion 
visible in the 12 CO J =2-1 maps bv lNakashima et al.l d2010t) at 
offsets of ~4 arcsec. However, it is noteworthy to remark that, 
specially in this case, we do not have enough observational data 
to further constrain the blobs geometry, since there is no indica- 
tion whatsoever of the spatial origin of the high-velocity emis- 
sion at higher J. Therefore, the location of the hotter and denser 
cylinder section in the aft (and not in the fore as a bow-shock) is 
purely arbitrary, and have no further consequences in the mod- 
eling. Given its slightly larger velocity, and its increased den- 
sity and temperature, it would be tempting to relate these jumps 
to the action of shocks, as a result of a front shock travelling 
through the nebula and catching up on the blobs, or a bow-shock 
in the fore section as the blobs expand against the PDR and 
molecular envelope. Unfortunately, the data at hand prevents fur- 
ther interpretation. We can only conclude on the existence of at 
least two components with different physical conditions. Further 
observations should provide independent insight on the proper- 
ties of these blobs. 

From our model we estimate a total mass of the blobs of 
0.1 M Q , under the aforementioned uncertainties. Therefore the 
total mass of the nebula, taking into account all the components, 
agrees well with previous determinations. 

Finally, at the adopted distance, the kinematical age of the 
blobs is in the range 360-630 yrs. Given that this value range 
is not model-dependent, this implies that the blobs are much 
younger than the rest of the molecular envelope, including the 
shocked, inner shell, whose kinematical age is a factor 3-4 larger. 
This suggests that the blobs were ejected approximately dur- 
ing the same epoch as (or immediately after than) the inner H 
ii n ebula, which s trongly emits at radio wavelengths (600 yrs, 
see M asson 1989). An alternative explanation, however, would 
be that the blobs are in fact older, but have been progressively 
accelerated by ram pressure by the faster-expanding inner H n 
nebula. 

4.3. Uncertainty of the model parameters 

Given the high-degree of interaction between the different model 
parameters, it is difficult to assess their range of uncertainty. 



Contrary to standard spatio-kinematical modeling, where the 
effect of changing a parameter is easily followed (e.g. a ve- 
locity near the lower limit of the uncertainty range corre- 
spo nds to an age-distance param eter close to the upper limit, 
see ISantander-Garcia et al. 1 120041 as an example), the effects 
of almost any parameter change in a radiative-transfer spatio- 
kinematical model are non-trivial to predict. 

Tables [2] and [3] show conservative estimations of the uncer- 
tainty range of every model parameter. The range for each pa- 
rameter has been obtained by fixing every other parameter of 
the model and allowing the given parameter to change until it no 
longer provided a fair fit to the data (i.e. until the intensity factors 
to be applied to the model profiles no longer lied within the 25% 
of the observed intensities, except in the case of the velocities, 
where the fairness of the fit was assessed by eye). Most of the 
uncertainty ranges of the model lie within a 10% of the best-fit 
values, except the temperature and density of the Ml shell, both 
of which reach a 20%; its abundance, which reaches a 40%; and 
the temperature, abundance and temperature of the blobs, some 
of which reach an uncertainty of ~40% of its best-fit value. 

These uncertainty ranges, however, should be considered as 
a simplistic approach to the real errors. In fact, if multiple pa- 
rameters were allowed to change, the uncertainty ranges would 
be somewhat larger: for example, an increase in temperature can 
be partially acommodated by a decrease in density, and vicev- 
ersa. Since, given the parameter interdependence, it would be 
impossible to compute the uncertainty ranges with as many de- 
grees of freedom as parameters in the model, we hereby pro- 
vide the following representative example. The temperature of 
the inner shell can be increased to 520 K if the density drops to 
1.3x10 s cm 4 , or the temperature decreased to 360 K and the 
density increased to 1.6xl0 5 cm 3 , and still provide a fair-fit to 
the data. This represents an increase in uncertainty for the tem- 
perature from 10% to 30%, and from 10% to 13% for the density. 

4.4. Other molecules 

The electronic level population is significantly more complicated 
for molecules more complex than CO, making the computation 
of their absorption and emission coefficients a daunting task. 
Hence, a detailed analysis including radiative-transfer solving 
of other molecular species such as H2O in NGC 7027 has been 
kept out of the focus of this study. However, it is interesting to do 
a qualitative analysis of the data in order to lay the foundations 
of future work. 
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Table 3. Best-fit model parameters for the high-velocity polar blobs of NGC 7027. The sizes correspond to a distance to the nebula 
of 1 kpc. The uncertainty range is included for every parameter except the distance of the structure to the central star, which is 
assumed a priori based on existing maps. 



Blob section 


Oil 


r f 


h* 


V* 


X( l2 CO) 


n 


T 




(10 16 cm) 


(10 16 cm) 


(10 1S cm) 


(kms 1 ) 




(cm 3 ) 


(K) 


Inner 


4.5 


3 +0 2 


1 65 +015 


40.0+! 


l^jxlO- 4 


1.5!^,xl0 5 


350 +25 


Outer 


7.65 


3 +0 3 


4 5 +03 


38.01! 


2.3^xl0- 4 


l^xlO 5 


50+ 25 

JU -1S 



Notes: The X( n CO) abundances are 50 times lower than their X( 12 CO) counterparts. 
' Distance from the star to the center of mass of the given cylindrical section. 

* Radius of the cylinder. 

* Height of the cylindrical section. 

* Velocity of the cylindrical section along the main axis of the nebula. 



H2O: The presence of the H2O li.o— lo,i an d li,i-Oo,o tran- 
sitions in the spectrum of NGC 7027 is striking, given that 
NGC 7027 is a C-rich PNe. In fact, to our knowledge, aside 
from this case, H2O has been de tected in only tw o C-rich PNe, 
CRL 618 dBuiarrabal et al.ll2012h and CRL 2688 dWesson et al.l 

[2oToh . 

The H2O profiles (see Fig. 5 in lBuiarrabal et al.l2012l) some- 
what resemble the high-excitation profiles of 12 CO (especially 
that of /=16-15). The relative weakness of these profiles, in 
comparison with those of CO at similar frequencies (i.e. simi- 
lar beam size), reveals that the nebula is optically thin at those 
transitions. The comparison of the H2O and CO profiles, in view 
of our model of the molecular envelope, seems to suggest that 
most of the water lies within the inner, faster shell, or even closer 
to the PDR, since the velocity separation between two peaks of 
the H2O li,o— lo,i an d 1 1 , 1 — Oo^o profiles are ~30.8 and ~31.9 
km s 'respectively. These values are slightly larger than the peak 
separation in every CO profile, including the highest /-line ob- 
served, 16-15, which is ~29.7 km s . 

It is known that water is alre ady abundant in the circumstel- 
lar envelopes around AGB stars dNeufeld et al.ll2.QlH) . However, 
the fact that water seems to be particularly abundant in the in- 
ner shells of NGC 7027 suggests an ongoing process of wa- 
ter production in this source, either due to the passage of a 
shock front or photo-induced chemistry. We consider photo- 
induction by the UV radiation from the central star to be the 
more likely cause, since while the shock is relatively weak (see 
section 4.2.1), most of the water would lie close to the ex- 
tremely hot (T e ff ~200,00 K) star. This explana tion is com- 
patible with the findings of iBujarrabal et al.l d2012l) in CRL 618 
and CRL 2688. CRL 618 shows a narrower water profile than 
those of CO, indicative of it being concentrated closer to the hot 
(TefT ~30,000 K) star, while there is no significant increase of 
the abundance in the shocked region. On the contrary, CRL 2688 
shows a shocked region but no water, while the central star is too 
cool (T e ff ~7,000 K) to emit a significant amount of UV radia- 
tion. 

C l& 0: The emission in C 18 Q is weak and optically thin (see 
Fig. 4 in lBuiarrabal et al.l2012l) . In view of our model, it appears 
to come primarily from the regions close to the equatorial plane. 
This could indicate that this ring is denser than the rest of the 
bipolar shell, or show matter condensations, in a very clumpy 
general structure. 

OH: The OH p rofile is very weak (see Fig. 5 in 
IBujarrabal et al.ll2012l) . Although OH is u sually associated with 
the re gion where H2O is dissociated (e.g. lGoldreich & Scovilld 
Il976l) . the shape of the profile is very different than those of wa- 



ter. In light of the present data, any further attempt at interpreting 
this profile is unclear. 

5. Conclusions 

The main results of this work are summarized in what follows: 

1 . We have developed a software code, sh apemol, which, used 
along with SHAPE (Stef fen et al.ll2.01 lb . implements spatio- 
kinematical modeling with accurate calculations of non-LTE 
line excitations and radiative transfer in molecular species. 
A more detailed description of shapemol will be given by 
Santander-Garcfa et al. in a forthcoming paper. 

2. Based on observations of a series of low and high-excitation 
12 CO and 13 CO transitions by Herschel/HIFI and the IRAM 
30-m radiotelescope, we have used shapemol and SHAPE to 
build a model of the molecular envelope of the young PNe 
NGC 7027. The geometry of the model is based in existing 
maps and data from CO and H2 (see section 4. 1). This model 
consists of four nested mildly-bipolar shells and a pair of 
high-velocity polar blobs (see Fig. [2j. Each shell is chara- 
terised by single values of the abundance, temperature, den- 
sity and velocity. Since the profiles contain little information 
on the spatial structure of NGC 7027, the geometry of our 
model should be considered as a rough approximation. 

3. The temperature of the shells ranges from 400 K for the 
innermost and 25 K for the outermost ones. The density 
also drops down from 1.5xl0 5 cm -3 in the innermost shell 
to 5xl0 3 cm -3 in the outermost one, while the 12 CO rel- 
ative abundance increases from 8xl0~ 5 to 6.7xl0 -4 . The 
12 CO/ 13 CO abundance ratio is 50, while the microturbu- 
lence, r5v, is 2 km s . 

4. The expansion pattern of the shells is not compatible with 
a Hubble-like flow. Instead, the 3-D velocity field of each 
shell is a combination of a constant, purely radial motion 
plus an additional constant component along the main axis 
of the nebula. The latter is only activated for distances to the 
equator larger than 1.5xl0 16 cm (at the adopted distance of 
1 kpc to the nebula). The velocity of each consecutive outer 
shell is slightly larger than the preceding inner one, with the 
following exception: 

5. The innermost shell is thinner than the rest, and is charac- 
terized by anomalous physical conditions: its expansion ve- 
locity is 33% higher than that in the adjacent intermediate 
shell. The abundance also increases a factor 1 .25 across the 
inner-to-intermediate shell discontinuity, while the density 
and temperature drop by factors of 3.75 and 4 respectively. 
We interpret these jumps as indicative of a passing shock 



s 



M. Santander-Garcfa et al.: Modeling the physical and excitation conditions of the molecular envelope of NGC 7027 



front. This shock could have important implications in the 
shaping of the nebula. 

6. Each of the two opposing polar blobs is assumed to be 
a small cylinder split in two sections with very different 
physical conditions. One of them is hot (350 K) and dense 
(l-5xl0 5 cm -3 ), has a 12 CO abundance of lxl(T 4 and a 
slightly larger velocity (40 km s _1 ) than the other one (38 
km s _1 ), which is much cooler (50 K), tenuous (1x10 s cm" 3 ) 
and shows a larger 12 CO abundance (2.3x10 4 ). The shape 
of the blobs, as well as the relative position of the cylinders 
(whether the hotter and denser is located at the fore or the aft 
of the blob) in the model is an arbitrary choice, given that the 
emission from the blobs is unresolved in the maps available 
in the literature and cannot therefore be determined from ex- 
isting data. Whether the blobs are affected by the same shock 
front as the inner shell or they form bow shocks in their ad- 
vance against the PDR and molecular envelope is therefore 
unclear. 

7. The computed molecular mass of the main nebula is 1.15 
M Q , while the mass of the blobs is 0. 1 M , adding to a to- 
tal molecular mass for NGC 7027 of roughly 1.3 M . This 
figure is compatible with those derived by previous works. 

8. The presence of H2O in the spectrum of NGC 7027, a C- 
rich nebula, is striking. From the profiles of the two detected 
transitions we conclude that the water is formed in a region 
close to the shocked inner shell, if not even slightly closer to 
the extremely hot central star. We cannot rule out the possi- 
bility of shocks as the leading mechanism in the formation 
of water. However, given the proximity to the extremely hot 
central star, we consider more likely that the formation is 
mainly photo-induced by UV radiation from the star. 
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Fig. 1. Resulting synthetic spectra (red, dotted-dashed line) and observations (black histogram) for the 12 CO and I3 CO transitions 
detected in NGC 7027. /factor refers to the intensity factor applied to the model to account for the uncertainties in the model's 
radiative transfer solving and in the observations' calibration. 
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Fig. 2. Left: Schematic representation of the model of the molecular envelope of NGC 7027. The shells and blobs are colored 
according to their densities, in a logarithmic scale. The axes denote the equatorial and polar directions. Right: Velocity field of 
the model, resulting from the combination of radial and polar components, both pointing outwards from the central star. The polar 
component is activated for distances to the equator greater than 1.5xl0 16 cm (the region where it is not active is indicated). The 
velocity of the polar blobs is purely along the polar direction. 
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